#version 420

layout(points) in;
layout(line_strip, max_vertices = 10) out;

uniform sampler3D TexDiagTensor;
uniform sampler3D TexOffDiagTensor;
uniform sampler3D TexMask;

uniform float     CurvatureThreshold;
uniform float     StopFA;
uniform float     StepSize;

uniform float     TotalLoop;
uniform vec3      KSizeScale;

/*
 |     boarder       |    inside      |      continue    |
0.0   0.166667   0.333333   0.5   0.666667   0.833333   1.0
*/
const float FinishTag = 0.8;

/*
 * \note:
 * tensors have to be stored in the following way:
 *
 * diag.x = Dxx     Diagonal elements are stored in one vector
 * diag.y = Dyy
 * diag.z = Dzz
 * offdiag.x = Dyz  Off diagonal elements are stored in another vector
 * offdiag.y = Dxz      where all components are stored in that location
 * offdiag.z = Dxy      that is indexed by the coordinate that is not in the tensor index
*/

float sqr(float a)
{
    return a * a;
}

vec3 getEigenvaluesCardano(vec3 diag, vec3 offdiag)
{
  const float M_SQRT3 = 1.73205080756887729352744634151;
  float de = offdiag.z * offdiag.x;
  float dd = sqr( offdiag.z );
  float ee = sqr( offdiag.x );
  float ff = sqr( offdiag.y );
  float m  = diag.x + diag.y + diag.z;
  float c1 = diag.x * diag.y + diag.x * diag.z + diag.y * diag.z
             - ( dd + ee + ff );
  float c0 = diag.z * dd + diag.x * ee + diag.y * ff - diag.x * diag.y * diag.z - 2. * offdiag.y * de;
  float p, sqrt_p, q, c, s, phi;
  p = sqr( m ) - 3. * c1;
  q = m * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
  sqrt_p = sqrt( abs( p ) );

  phi = 27. * ( 0.25 * sqr( c1 ) * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
  phi = ( 1. / 3. ) * atan( sqrt( abs( phi ) ), q );

  c = sqrt_p * cos( phi );
  s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );

  vec3 w;
  w[2] = ( 1. / 3. ) * ( m - c );
  w[1] = w[2] + s;
  w[0] = w[2] + c;
  w[2] -= s;
  return w;
}

vec3 getEigenvector(vec3 ABC /*diag without eigenalue i*/, vec3 offdiag)
{
  vec3 vec;
  vec.x = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.x - ABC.z * offdiag.z );
  vec.y = ( offdiag.y * offdiag.x - ABC.z * offdiag.z ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );
  vec.z = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );

  return normalize( vec );
}

float getFA( vec3 diag, vec3 offdiag )
{
  float cross = dot( offdiag, offdiag );
  float j2 = diag.x * diag.y + diag.y * diag.z + diag.z * diag.x - cross;
  float j4 = dot( diag, diag ) + 2. * cross;
  return sqrt( ( j4 - j2 ) / j4 );
}

void getPrimaryEigenVectorAndFA(vec3 pos, out vec3 eVector, out float FA)
{
    // this is in model coordinate space
    vec3 texCoord = pos * KSizeScale;
    vec3 diagnal = texture(TexDiagTensor, texCoord).xyz;
    vec3 offDiag = texture(TexOffDiagTensor, texCoord).xyz;

    vec3 eigenValue = getEigenvaluesCardano(diagnal, offDiag);

    eVector = getEigenvector(diagnal-eigenValue.x, offDiag);
    FA = getFA(diagnal, offDiag);
}

struct Vertex
{
    vec4 pos;	// the unit should be mm
    vec3 direction;
    float FA;
};

void main()
{
    bool isTrackingEnd = false;
    Vertex pcurt = Vertex(vec4(0.0), vec3(0.0), 0.0);   // current point
    Vertex pnext = Vertex(vec4(0.0), vec3(0.0), 0.0);   // next point

    pcurt.pos = gl_in[0].gl_Position;

    getPrimaryEigenVectorAndFA(pcurt.pos.xyz, pcurt.direction, pcurt.FA);
    int loop = 0;
    while((!isTrackingEnd) && (TotalLoop > loop))
    {
        ++loop;
        
        // make a small step at direction of the current point primary eigenvector
        pnext.pos = vec4(pcurt.pos.xyz + 
            sign(pcurt.pos.w) * StepSize * pcurt.direction, pcurt.pos.w);

        //1. test if current point inside the mask region
        //float mask = texture(TexMask, currentPosition).x;

        
        if(!isTrackingEnd)
        {
            getPrimaryEigenVectorAndFA(pnext.pos.xyz, pnext.direction, pnext.FA);
            if(StopFA > pnext.FA)
            {
                isTrackingEnd = true;
            }
        }
        
        // calculate the curvature between the current point primary eigenvector and
        // the next point primary eigenvector, and test the curvature if it exceed the
        // the CurvatureThreshold threshold.
        if(!isTrackingEnd)
        {
            float curvature = abs(dot(pnext.direction, pcurt.direction));
            if(CurvatureThreshold > curvature) // the angle is to larger
            {
                isTrackingEnd = true;
            }
        }
        
        if(isTrackingEnd)
        {
            pcurt.pos.w += sign(pcurt.pos.w) * FinishTag;
        }
        
        gl_Position = pcurt.pos;
        EmitVertex();
        
        pcurt = pnext;
    }
    EndPrimitive();
}	
    //bool bTrackingEnd = false;
//
    //int Iteration = 0;
    //float cosfMaxBendAngle = cos(fMaxBendAngle);
//
    //// vec4 v4CurCoord = gl_PositionIn[0];
    //vec4 v4CurCoord = vec4(0.0);
    //vec4 v4NextPoint = vec4(0.0);
//
    //vec3 v3Diagonal = vec3(0.0);
    //vec3 v3OffDiag  = vec3(0.0);
//
    //vec3 v3EigenValues = vec3(0.0);
//
    //vec3 v3CurEigenVector = vec3(0.0);
    //vec3 v3NextEigenVector = vec3(0.0);
//
    //v4CurCoord = gl_in[0].gl_Position;
    //v4CurCoord.xyz += vec3(0.5,0.5,0.5);
    //v4CurCoord.xyz /= v3VolumeSize;
//
    //v3Diagonal = texture(TensorDiagonal,v4CurCoord.xyz).xyz;
    //v3OffDiag  = texture(TensorOffDiagonal,v4CurCoord.xyz).xyz;
    //v3EigenValues = getEigenvaluesCardano(v3Diagonal,v3OffDiag);
//
    //v3CurEigenVector = getEigenvector(v3Diagonal-v3EigenValues.xxx,v3OffDiag);
//
    //while(!bTrackingEnd && Iteration < iterationTimes)
    //{
        //if(all(lessThan(abs(v3CurEigenVector),vec3(0.001,0.001,0.001))))
        ////if(any(lessThan(abs(v3CurEigenVector),vec3(0.001,0.001,0.001))))
        //{
            //bTrackingEnd = true;
        //}
//
        //if(!bTrackingEnd)
        //{
            //// step to the next adjecent point
            //v4NextPoint = vec4(
                //v4CurCoord.xyz + sign(v4CurCoord.w) * fSampleStride / v3VolumeSize * v3CurEigenVector,
                //v4CurCoord.w);
            //// bvec3 less = lessThan(v4NextPoint.xyz, vec3(0.5,0.5,0.5)/ VolumeSize);
            //// bvec3 big  = greaterThan(v4NextPoint.xyz, vec3(1.0,1.0,1.0) -  vec3(0.5,0.5,0.5) / VolumeSize);
//
            //bvec3 less = lessThan(v4NextPoint.xyz, vec3(0.0));
            //bvec3 big =  greaterThan(v4NextPoint.xyz, vec3(1.0));
//
            //// if the next point exceed the texture border
            //if(any(less)||any(big))
            //{
                //bTrackingEnd = true;
            //}
        //}
        //
        //if(!bTrackingEnd)
        //{
            ////FiberColor = vec4(ColorEncoding1(v3EigenValues),0.5);
            //
            //// calculate the next point eigenvector
            //v3Diagonal = texture(TensorDiagonal,v4NextPoint.xyz).xyz;
            //v3OffDiag  = texture(TensorOffDiagonal,v4NextPoint.xyz).xyz;
            //v3EigenValues = getEigenvaluesCardano(v3Diagonal,v3OffDiag);
//
            //v3NextEigenVector = getEigenvector(v3Diagonal-v3EigenValues.xxx,v3OffDiag);
//
            //float cosAngle = dot(v3NextEigenVector, v3CurEigenVector);
//
            //if((length(v3Diagonal) == 0.0) || (cosfMaxBendAngle > abs(cosAngle))) //todo
            //{
                //bTrackingEnd = true;
            //}
        //}
//
        //outPoint = vec4(v4CurCoord.xyz * v3VolumeSize, v4CurCoord.w);
//
        //if(bTrackingEnd)
        //{
            //outPoint.w = trunc(v4CurCoord.w) + fBorderTag * sign(v4CurCoord.w);
        //}
        //else
        //{
            //// is last iteration?
            //if(iterationTimes - Iteration < 1.1)
            //{
                //outPoint.w =trunc(v4CurCoord.w) + fContinueTag * sign(v4CurCoord.w);
            //}
            //else
            //{
                //outPoint.w = trunc(v4CurCoord.w) + fInsideTag * sign(v4CurCoord.w);
            //}
        //}
        //EmitVertex();
        //EndPrimitive();
//
        //Iteration++;
        //v4CurCoord = v4NextPoint;
        //v3CurEigenVector = v3NextEigenVector;
    //}
//}
